Paso 5 (1)- Trayectorias de hospitalización y mortalidad con foco en condiciones vinculadas a trastornos de salud mental y consumo de sustancias posterior a un primer ingreso por alguno de estos trastornos, en usuarios/as jóvenes y adultos emergentes de población general y pertenecientes a pueblos originarios, 2018-2021, Chile
Representar las mejores opciones de agrupamiento, junto con su relación con otras variables. En esta etapa, analizamos la asociación con covariables de la solución con mejores índices de calidad para la resolución trimestral.
Autor/a
Andrés González Santa Cruz
Fecha de publicación
3 de abr, 2025
Configurar
Código
# remover objetos y memoria utilizadarm(list=ls());gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 599277 32.1 1288917 68.9 690014 36.9
Vcells 1144548 8.8 8388608 64.0 1879001 14.4
#elegir repositorioif(Sys.info()["sysname"]=="Windows"){options(repos =c(CRAN ="https://cran.dcc.uchile.cl/"))}options(install.packages.check.source ="yes") # Chequea la fuente de los paquetes#borrar caché#system("fc-cache -f -v")if(!require(pacman)){install.packages("pacman");require(pacman)}pacman::p_unlock(lib.loc =.libPaths()) #para no tener problemas reinstalando paquetesif(Sys.info()["sysname"]=="Windows"){if (getRversion() !="4.4.1") { stop("Requiere versión de R 4.4.1. Actual: ", getRversion()) }}cat("quarto version: "); system("quarto --version")
quarto version:
[1] 0
Código
if(!require(job)){install.packages("job");require(job)}if(!require(kableExtra)){install.packages("kableExtra");require(kableExtra)}if(!require(tidyverse)){install.packages("tidyverse");require(tidyverse)}if(!require(cluster)){install.packages("cluster"); require(cluster)}if(!require(WeightedCluster)){install.packages("WeightedCluster"); require(WeightedCluster)}if(!require(devtools)){install.packages("devtools"); require(devtools)}if(!require(TraMineR)){install.packages("TraMineR"); require(TraMineR)}if(!require(TraMineRextras)){install.packages("TraMineRextras"); require(TraMineRextras)}if(!require(NbClust)){install.packages("NbClust"); require(NbClust)}if(!require(haven)){install.packages("haven"); require(haven)}if(!require(ggseqplot)){install.packages("ggseqplot"); require(ggseqplot)}if(!require(grid)){install.packages("grid"); require(grid)}if(!require(gridExtra)){install.packages("gridExtra"); require(gridExtra)}if(!require(Tmisc)){install.packages("Tmisc"); require(Tmisc)}if(!require(factoextra)){install.packages("factoextra"); require(factoextra)}if(!require(stargazer)){install.packages("stargazer"); require(stargazer)}if(!require(gtsummary)){install.packages("gtsummary"); require(gtsummary)}if(!require(lmtest)){install.packages("lmtest"); require(lmtest)}if(!require(emmeans)){install.packages("emmeans"); require(emmeans)}if(!require(fpp2)){install.packages("fpp2"); require(fpp2)}if(!require(purrr)){install.packages("purrr"); require(purrr)}if(!require(forecast)){install.packages("forecast"); require(forecast)}if(!require(magrittr)){install.packages("magrittr"); require(magrittr)}if(!require(foreach)){install.packages("foreach"); require(foreach)}if(!require(doParallel)){install.packages("doParallel"); require(doParallel)}if(!require(progressr)){install.packages("progressr"); require(progressr)}if(!require(chisq.posthoc.test)){devtools::install_github("ebbertd/chisq.posthoc.test")}if(!require(rstatix)){install.packages("rstatix"); require(rstatix)}if(!require(rio)){install.packages("rio"); require(rio)}if(!require(cowplot)){install.packages("cowplot"); require(cowplot)}if(!require(DiagrammeR)){install.packages("DiagrammeR"); require(DiagrammeR)}if(!require(DiagrammeRsvg)){install.packages("DiagrammeRsvg"); require(DiagrammeRsvg)}if(!require(rsvg)){install.packages("rsvg"); require(rsvg)}if(!require(survminer)){install.packages("survminer"); require(survminer)}if(!require(epitools)){install.packages("epitools"); require(epitools)}seq_mean_t_dos_grupos <-function(bd =NULL, group1, group2) {# Agrupar por ambas variables resultados <-by(bd, list(group1, group2), seqmeant)# Obtener todas las combinaciones posibles de los grupos combinaciones <-expand.grid(group1 =unique(group1), group2 =unique(group2), stringsAsFactors =FALSE)# Extraer los resultados y asociarlos con las combinaciones resultados_df <-do.call(rbind, lapply(seq_along(resultados), function(i) { group_name1 <-attr(resultados, "dimnames")[[1]][i] group_name2 <-attr(resultados, "dimnames")[[2]][i]data.frame(factor_inclusivo_1 = group_name1, factor_inclusivo_2 = group_name2, Mean = resultados[[i]]) }))# Unir los resultados con las combinaciones para rellenar los valores faltantes final_df <-merge(combinaciones, resultados_df, by.x =c("group1", "group2"), by.y =c("factor_inclusivo_1", "factor_inclusivo_2"), all.x =TRUE)return(final_df)}multinom_pivot_wider <-function(x) {# check inputs match expectatations# create tibble of results df <- tibble::tibble(outcome_level =unique(x$table_body$groupname_col)) df$tbl <- purrr::map( df$outcome_level,function(lvl) { gtsummary::modify_table_body( x, ~dplyr::filter(.x, .data$groupname_col %in% lvl) |> dplyr::ungroup() |> dplyr::select(-.data$groupname_col) ) } )tbl_merge(df$tbl, tab_spanner =paste0("**", df$outcome_level, "**"))}best_subset_multinom <-function(y, x.vars, data) {# y Nombre de la variable dependiente (cadena de texto)# x.vars Vector de nombres de predictores (caracter)# data Dataframe con los datos de entrenamiento# Cargar las librerías necesariasrequire(dplyr)require(purrr)require(tidyr)require(nnet)require(MASS)# Generar todas las combinaciones posibles de predictores predictors_list <-lapply(1:length(x.vars), function(i) {combn(x.vars, i, simplify =FALSE) }) |>unlist(recursive =FALSE)# Inicializar una lista para almacenar los resultados results <-list()# Iterar sobre cada combinación de predictoresfor (i inseq_along(predictors_list)) { predictors <- predictors_list[[i]] formula <-as.formula(paste(y, "~", paste(predictors, collapse ="+")))# Ajustar el modelo multinomial model <-tryCatch( nnet::multinom(formula, data = data, trace =FALSE),error =function(e) NULL )# Si el modelo se ajustó correctamente, almacenar los resultadosif (!is.null(model)) {# Extraer el AIC del modelo aic <-AIC(model)# Almacenar la información en una lista results[[length(results) +1]] <-list(predictors = predictors,model = model,AIC = aic ) } }# Convertir la lista de resultados en un dataframe results_df <- results |> purrr::map_df(function(res) {data.frame(predictors =paste(res$predictors, collapse ="+"),AIC = res$AIC,stringsAsFactors =FALSE ) })# Ordenar los modelos por AIC de menor a mayor results_df <- results_df |>arrange(AIC)return(results_df)}best_subset_multinom_interactions <-function(y, x.vars, data) {# y Nombre de la variable dependiente (cadena de texto)# x.vars Vector de nombres de predictores (caracter)# data Dataframe con los datos de entrenamiento# Cargar las librerías necesariasrequire(dplyr)require(purrr)require(tidyr)require(nnet)require(MASS)# Generar todas las combinaciones posibles de predictores (efectos principales) main_effects_list <-lapply(1:length(x.vars), function(i) {combn(x.vars, i, simplify =FALSE) }) |>unlist(recursive =FALSE)# Inicializar una lista para almacenar los resultados results <-list()# Iterar sobre cada combinación de efectos principalesfor (main_effects in main_effects_list) {# Generar términos de interacción de hasta 3 variables interaction_terms <-list()# Para interacciones de 2 variablesif (length(main_effects) >=2) { interaction_terms_2way <-combn(main_effects, 2, function(x) paste(x, collapse =":")) interaction_terms <-c(interaction_terms, interaction_terms_2way) }# Para interacciones de 3 variablesif (length(main_effects) >=3) { interaction_terms_3way <-combn(main_effects, 3, function(x) paste(x, collapse =":")) interaction_terms <-c(interaction_terms, interaction_terms_3way) }# Combinar efectos principales e interacciones all_terms <-c(main_effects, interaction_terms)# Generar todas las combinaciones posibles de términos (incluyendo interacciones)# Solo se incluyen interacciones si sus efectos principales están presentes term_combinations <-list()# Obtener todos los subconjuntos de efectos principales main_effects_subsets <-lapply(1:length(main_effects), function(i) {combn(main_effects, i, simplify =FALSE) }) |>unlist(recursive =FALSE)# Para cada subconjunto de efectos principalesfor (me in main_effects_subsets) {# Iniciar con los efectos principales terms <- me# Incluir interacciones solo si todos sus efectos principales están incluidos possible_interactions <- interaction_terms[sapply(interaction_terms, function(x) { vars_in_interaction <-unlist(strsplit(x, ":"))all(vars_in_interaction %in% me) }) ]# Generar todas las combinaciones de interacciones para incluir interaction_subsets <-list(NULL)if (length(possible_interactions) >0) { interaction_subsets <-lapply(1:length(possible_interactions), function(i) {combn(possible_interactions, i, simplify =FALSE) }) |>unlist(recursive =FALSE) }# Para cada combinación de interacciones, crear el conjunto completo de términosfor (ints in interaction_subsets) {if (is.null(ints)) { full_terms <- terms } else { full_terms <-c(terms, ints) }# Añadir a la lista de combinaciones de términos term_combinations <-append(term_combinations, list(full_terms)) } }# Ajustar modelos para cada combinación de términosfor (terms in term_combinations) { formula <-as.formula(paste(y, "~", paste(terms, collapse ="+")))# Ajustar el modelo multinomial model <-tryCatch( nnet::multinom(formula, data = data, trace =FALSE),error =function(e) NULL,warning =function(w) NULL )# Si el modelo se ajustó correctamente, almacenar los resultadosif (!is.null(model)) {# Extraer el BIC del modelo bic <-BIC(model)# Almacenar la información en la lista de resultados results[[length(results) +1]] <-list(predictors =paste(terms, collapse =" + "),model = model,BIC = bic ) } } }# Convertir la lista de resultados en un dataframe results_df <- results |> purrr::map_df(function(res) {data.frame(predictors = res$predictors,BIC = res$BIC,stringsAsFactors =FALSE ) })# Ordenar los modelos por BIC de menor a mayor results_df <- results_df |>arrange(BIC)return(results_df)}best_subset_multinom_interactions_parallel <-function(y, x.vars, data) {# y Nombre de la variable dependiente (cadena de texto)# x.vars Vector de nombres de predictores (caracter)# data Dataframe con los datos de entrenamiento# Cargar las librerías necesarias dentro de la funciónrequire(dplyr)require(purrr)require(tidyr)require(nnet)require(MASS)require(foreach)require(doParallel)require(progressr)# Iniciar los gestores de progresohandlers(global =TRUE)handlers("txt")# Generar todas las combinaciones posibles de predictores (efectos principales) main_effects_list <-lapply(1:length(x.vars), function(i) {combn(x.vars, i, simplify =FALSE) }) |>unlist(recursive =FALSE)# Inicializar una lista para almacenar las fórmulas de los modelos formulas_list <-list()# Generar todas las fórmulas posibles con interacciones hasta de 3 variablesfor (main_effects in main_effects_list) {# Generar términos de interacción de hasta 3 variables interaction_terms <-character(0) # Aseguramos que es un vector de caracteres# Para interacciones de 2 variablesif (length(main_effects) >=2) { interaction_terms_2way <-combn(main_effects, 2, function(x) paste(x, collapse =":"), simplify =TRUE) interaction_terms <-c(interaction_terms, interaction_terms_2way) }# Para interacciones de 3 variablesif (length(main_effects) >=3) { interaction_terms_3way <-combn(main_effects, 3, function(x) paste(x, collapse =":"), simplify =TRUE) interaction_terms <-c(interaction_terms, interaction_terms_3way) }# Generar todas las combinaciones posibles de efectos principales main_effects_subsets <-lapply(1:length(main_effects), function(i) {combn(main_effects, i, simplify =FALSE) }) |>unlist(recursive =FALSE)# Para cada subconjunto de efectos principalesfor (me in main_effects_subsets) {# Iniciar con los efectos principales terms <- me# Identificar interacciones cuyos efectos principales están en 'me'if (length(interaction_terms) >0) { possible_interactions <- interaction_terms[vapply(interaction_terms, function(x) { vars_in_interaction <-unlist(strsplit(x, ":"))all(vars_in_interaction %in% me) }, FUN.VALUE =logical(1)) ] } else { possible_interactions <-character(0) }# Generar todas las combinaciones posibles de estas interacciones interaction_subsets <-list(character(0)) # Incluir el caso sin interaccionesif (length(possible_interactions) >0) { interaction_combinations <-lapply(1:length(possible_interactions), function(i) {combn(possible_interactions, i, simplify =FALSE) }) |>unlist(recursive =FALSE) interaction_subsets <-c(interaction_subsets, interaction_combinations) }# Para cada combinación de interaccionesfor (ints in interaction_subsets) { full_terms <-c(terms, ints)# Crear la fórmula del modelo y almacenarla formula_str <-paste(y, "~", paste(full_terms, collapse ="+")) formulas_list <-append(formulas_list, list(formula_str)) } } }# Eliminar posibles duplicados de fórmulas formulas_list <-unique(formulas_list)# Total de modelos a ajustar total_models <-length(formulas_list)# Iniciar el progreso p <-progressor(steps = total_models)# Ajustar los modelos en paralelo usando foreach results_list <-foreach(i =1:total_models, .packages =c("nnet", "MASS"), .combine ='rbind') %dopar% { formula_str <- formulas_list[[i]] formula <-as.formula(formula_str)# Ajustar el modelo model <-tryCatch( nnet::multinom(formula, data = data, trace =FALSE),error =function(e) NULL,warning =function(w) NULL )# Actualizar el progresop(sprintf("Ajustando modelo %d de %d", i, total_models))# Si el modelo se ajustó correctamente, almacenar los resultadosif (!is.null(model)) { bic <-BIC(model)data.frame(predictors = formula_str,BIC = bic,stringsAsFactors =FALSE ) } else {NULL } }# Convertir los resultados a dataframe y ordenar por BIC results_df <-as.data.frame(results_list) results_df <- results_df |>arrange(BIC)return(results_df)}num_cores <- parallel::detectCores() -1cl <-makeCluster(num_cores)registerDoParallel(cl)#pacman job kableExtra tidyverse cluster WeightedCluster devtools TraMineR TraMineRextras NbClust haven ggseqplot gridExtra Tmisc factoextra reticulate withr rmarkdown quartooptions(knitr.kable.NA ='')#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:##:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#knitr::knit_hooks$set(time_it =local({ now <-NULLfunction(before, options) {if (before) {# record the current time before each chunk now <<-Sys.time() } else {# calculate the time difference after a chunk res <-ifelse(difftime(Sys.time(), now)>(60^2),difftime(Sys.time(), now)/(60^2),difftime(Sys.time(), now)/(60^1))# return a character string to show the time x<-ifelse(difftime(Sys.time(), now)>(60^2),paste("Tiempo que demora esta sección:", round(res,1), "horas"),paste("Tiempo que demora esta sección:", round(res,1), "minutos"))paste('<div class="message">', gsub('##', '\n', x),'</div>', sep ='\n') } }}))knitr::opts_chunk$set(time_it =TRUE)#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:format_cells <-function(df, rows ,cols, value =c("italics", "bold", "strikethrough")){# select the correct markup# one * for italics, two ** for bold map <-setNames(c("*", "**", "~~"), c("italics", "bold", "strikethrough")) markup <- map[value] for (r in rows){for(c in cols){# Make sure values are not factors df[[c]] <-as.character( df[[c]])# Update formatting df[r, c] <-ifelse(nchar(df[r, c])==0,"",paste0(markup, gsub(" ", "", df[r, c]), markup)) } }return(df)}#To produce line breaks in messages and warningsknitr::knit_hooks$set(error =function(x, options) {paste('\n\n<div class="alert alert-danger" style="font-size: small !important;">',gsub('##', '\n', gsub('^##\ Error', '**Error**', x)),'</div>', sep ='\n') },warning =function(x, options) {paste('\n\n<div class="alert alert-warning" style="font-size: small !important;">',gsub('##', '\n', gsub('^##\ Warning:', '**Warning**', x)),'</div>', sep ='\n') },message =function(x, options) {paste('<div class="message" style="font-size: small !important;">',gsub('##', '\n', x),'</div>', sep ='\n') })#_#_#_#_#_#_#_#_#_#_#_#_#_invisible("Function to format CreateTableOne into a database")as.data.frame.TableOne <-function(x, ...) {capture.output(print(x,showAllLevels =TRUE, varLabels = T,...) -> x) y <-as.data.frame(x) y$characteristic <- dplyr::na_if(rownames(x), "") y <- y |>fill(characteristic, .direction ="down") |> dplyr::select(characteristic, everything())rownames(y) <-NULL y}#_#_#_#_#_#_#_#_#_#_#_#_#_# Austin, P. C. (2009). The Relative Ability of Different Propensity # Score Methods to Balance Measured Covariates Between # Treated and Untreated Subjects in Observational Studies. Medical # Decision Making. https://doi.org/10.1177/0272989X09341755smd_bin <-function(x,y){ z <- x*(1-x) t <- y*(1-y) k <-sum(z,t) l <- k/2return((x-y)/sqrt(l))}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:format_table_vec <-function(tbl, digits =1) { counts <-as.numeric(tbl) percentages <-prop.table(tbl) *100 formatted <-sapply(seq_along(counts), function(i) { p_val <- percentages[i]# Si el porcentaje es prácticamente entero, formatea sin decimalesif (abs(p_val -round(p_val)) < .Machine$double.eps^0.5) { p_str <-sprintf("%.0f", p_val) } else { p_str <-sprintf(paste0("%.", digits, "f"), p_val) }paste0(counts[i], " (", p_str, ")") }) formatted}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:if(.Platform$OS.type =="windows") withAutoprint({memory.size()memory.size(TRUE)memory.limit()})
> memory.size()
[1] Inf
> memory.size(TRUE)
[1] Inf
> memory.limit()
[1] Inf
Código
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:func_tab_range_clus<-function(range_clus){rbind.data.frame(lapply(list(as.vector(rev(sort(table(range_clus$clustering$cluster2)))),as.vector(rev(sort(table(range_clus$clustering$cluster3)))),as.vector(rev(sort(table(range_clus$clustering$cluster4)))),as.vector(rev(sort(table(range_clus$clustering$cluster5)))),as.vector(rev(sort(table(range_clus$clustering$cluster6)))),as.vector(rev(sort(table(range_clus$clustering$cluster7)))),as.vector(rev(sort(table(range_clus$clustering$cluster8)))),as.vector(rev(sort(table(range_clus$clustering$cluster9)))),as.vector(rev(sort(table(range_clus$clustering$cluster10)))),as.vector(rev(sort(table(range_clus$clustering$cluster11)))),as.vector(rev(sort(table(range_clus$clustering$cluster12)))),as.vector(rev(sort(table(range_clus$clustering$cluster13)))),as.vector(rev(sort(table(range_clus$clustering$cluster14)))),as.vector(rev(sort(table(range_clus$clustering$cluster15)))) ),function(x) { length_out <-max(sapply(list(as.vector(rev(sort(table(range_clus$clustering$cluster2)))),as.vector(rev(sort(table(range_clus$clustering$cluster3)))),as.vector(rev(sort(table(range_clus$clustering$cluster4)))),as.vector(rev(sort(table(range_clus$clustering$cluster5)))),as.vector(rev(sort(table(range_clus$clustering$cluster6)))),as.vector(rev(sort(table(range_clus$clustering$cluster7)))),as.vector(rev(sort(table(range_clus$clustering$cluster8)))),as.vector(rev(sort(table(range_clus$clustering$cluster9)))),as.vector(rev(sort(table(range_clus$clustering$cluster10)))),as.vector(rev(sort(table(range_clus$clustering$cluster11)))),as.vector(rev(sort(table(range_clus$clustering$cluster12)))),as.vector(rev(sort(table(range_clus$clustering$cluster13)))),as.vector(rev(sort(table(range_clus$clustering$cluster14)))),as.vector(rev(sort(table(range_clus$clustering$cluster15)))) ), length))c(x, rep(NA, length_out -length(x))) } ))|>t() |>data.frame()|>`rownames<-`(NULL)}frobenius_norm <-function(matrix1, matrix2) {if (!all(dim(matrix1) ==dim(matrix2))) {stop("Matrices must have the same dimensions") }# Replace NA values with 0 (or any other desired default) matrix1[is.na(matrix1)] <-0 matrix2[is.na(matrix2)] <-0# Calculate the residuals residuals <- matrix1 - matrix2# Frobenius norm frobenius <-sqrt(sum(residuals^2))return(frobenius)}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:confcqi2 <-function(nullstat, quant, n){ alpha <- (1-quant)/2#calpha <- alpha+(alpha-1)/n#print(c(calpha, alpha))#minmax <- quantile(nullstat, c(calpha, 1-calpha)) minmax <-quantile(nullstat, c(alpha, 1-alpha))return(minmax)}normstatcqi2 <-function(bcq, stat, norm=TRUE){ origstat <- bcq$clustrange$stats[, stat] nullstat <- bcq$stats[[stat]]#normstat <- rbind(nullstat, origstat)if(norm){for(i inseq_along(origstat)){ mx <-mean(nullstat[, i]) sdx <-sd(nullstat[, i]) nullstat[ , i] <- (nullstat[, i]-mx)/sdx origstat[i] <- (origstat[i]-mx)/sdx } } alldatamax <-apply(nullstat, 1, max)#as.vector(xx) sumcqi <-list(origstat=origstat, nullstat=nullstat, alldatamax=alldatamax)return(sumcqi)}print.seqnullcqi.powder <-function(x, norm =FALSE, quant =0.95, digits =2, append =FALSE, ...) {cat("Parametric bootstrap cluster analysis validation\n")cat("Sequence analysis null model:", deparse(x$nullmodel), "\n")cat("Number of bootstraps:", x$R, "\n")cat("Clustering method:", ifelse(x$kmedoid, "PAM/K-Medoid", paste0("hclust with ", x$hclust.method)), "\n")cat("Seqdist arguments:", deparse(x$seqdist.args), "\n\n\n") alls <-as.data.frame(x$clustrange$stats) quants <-rep("", ncol(alls))names(quants) <-colnames(alls)for (ss incolnames(alls)) { sumcqi <-normstatcqi2(x, stat = ss, norm = norm) alls[, ss] <-as.character(round(sumcqi$origstat, digits = digits)) borne <-as.character(round(confcqi2(sumcqi$alldatamax, quant, x$R), digits = digits)) quants[ss] <-paste0("[", borne[1], "; ", borne[2], "]") } results_tibble <- tibble::as_tibble(rbind(alls, rep("", length(quants)), quants))# Print a summary to the console for immediate feedbackrownames(results_tibble) <-c(rownames(x$clustrange$stats), "", paste("Null Max-T", quant, "interval")) results_df <-as.data.frame(results_tibble)print(results_tibble, ...)return(list(results_tibble= results_tibble, results_df= results_df ))}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:##:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:# Función para aplicar la prueba de Fisher a todas las combinaciones de filas usando todas las columnasfisher_posthoc_all_cols <-function(contingency_table) {# Obtener combinaciones de filas (pares) row_pairs <-combn(rownames(contingency_table), 2, simplify =FALSE)# Aplicar la prueba de Fisher a cada par de filas usando todas las columnas al mismo tiempo results <-map_dfr(row_pairs, function(pair) {# Crear tabla de 2xN para el par de filas en todas las columnas sub_table <- contingency_table[pair, , drop =FALSE]# Aplicar el test de Fisher test_result <-fisher.test(sub_table, simulate.p.value=T,B=1e4)# Devolver los resultados en un data frametibble(Row1 = pair[1],Row2 = pair[2],p.value = test_result$p.value ) })# Ajustar p-valores usando el método de Holm results <- results |>mutate(p.adjusted =p.adjust(p.value, method ="holm"))return(results)}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:save_base_plot_as_grob <-function(plot_expr, res=300, width =1600, height=1200) {# Crea un archivo temporal con extensión .png filename <-tempfile(fileext =".png")# Guarda el gráfico en alta resolución en el archivo temporalpng(filename, width = width, height = height, res = res)replayPlot(plot_expr) # Reproduce el gráfico grabadodev.off() # Cierra el dispositivo gráfico# Convierte el archivo PNG en un objeto gráfico (grob) grob <- grid::rasterGrob(png::readPNG(filename), interpolate =TRUE)return(grob) # Devuelve el grob}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:chisq_cramerv<-function(contingency_table){ chisq_test <-chisq.test(contingency_table) cramers_v <-sqrt(chisq_test$statistic / (sum(contingency_table) * (min(dim(contingency_table)) -1)))list(chisq_statistic=sprintf("%1.2f", chisq_test$statistic), chisq_df= chisq_test$parameter, chisq_p_value =ifelse(chisq_test$p.value<.001, "<0.001", sprintf("%1.4f", chisq_test$p.value)), cramers_v =sprintf("%1.2f", cramers_v))}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#oneway_anova_effect_size <-function(values, group) {# Perform one-way ANOVA anova_result <-aov(values ~ group)# Summarize ANOVA results anova_summary <-summary(anova_result)# Extract sums of squares ss_between <- anova_summary[[1]]$"Sum Sq"[1] ss_total <-sum(anova_summary[[1]]$"Sum Sq")# Calculate eta-squared eta_squared <- ss_between / ss_total# Return ANOVA summary and effect sizelist(anova_summary = anova_summary,eta_squared = eta_squared )}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#duración en cada estadoseq_mean_t<-function(bd=NULL, group= group){ resultados<-by(bd, group, seqmeant)do.call(rbind, lapply(names(resultados), function(name) {data.frame(factor_inclusivo = name, resultados[[name]]) }))}seqtrate_t<-function(bd=NULL, group= group){# Utilizar la función 'by' para calcular las tasas agrupadas por 'glosa_sexo' resultados <-by(bd, group, seqtrate)# Convertir los resultados en un data frame en formato largo resultados_long <-do.call(rbind, lapply(names(resultados), function(sexo) { df <-as.data.frame(resultados[[sexo]]) df$from <-rownames(df) df$glosa_sexo <- sexo df }))# Usar tidyr para convertir a formato largolibrary(tidyr) resultados_long <-pivot_longer(resultados_long, cols =-c(from, glosa_sexo), names_to ="to", values_to ="rate")# Mostrar el data frame finalprint(resultados_long)}seqcount_t<-function(bd=NULL, group= group){# Utilizar la función 'by' para calcular las tasas agrupadas por 'glosa_sexo' resultados <-by(bd, group, function(x) seqtrate(x, count =TRUE))# Convertir los resultados en un data frame en formato largo resultados_long <-do.call(rbind, lapply(names(resultados), function(sexo) { df <-as.data.frame(resultados[[sexo]]) df$from <-rownames(df) df$glosa_sexo <- sexo df }))# Usar tidyr para convertir a formato largolibrary(tidyr) resultados_long <-pivot_longer(resultados_long, cols =-c(from, glosa_sexo), names_to ="to", values_to ="count")# Mostrar el data frame finalprint(resultados_long)}#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
Descartamos objetos que no serán de interés para el análisis
Código
# Define el vector con los nombres de los objetos que quieres conservarmis_objetos <-c("as.data.frame.TableOne", "avs_por_cluster_month", "avs_por_cluster_quarter","best_subset_multinom", "best_subset_multinom_interactions", "best_subset_multinom_interactions_parallel", "categories_hac_om2_m", "categories_pam_om4_q", "chisq_cramerv", "cl", "costmatrix_month", "costmatrix_quarter", "data_long_establecimiento_2024_std", "df_filled2", "dist_month_lcs", "dist_month_om", "dist_quarter_lcs", "dist_quarter_om", "dt_ing_calendar_month_t_desde_primera_adm_dedup", "dt_ing_calendar_quarter_t_desde_primera_adm_dedup", "fisher_posthoc_all_cols", "folder_path", "format_cells", "format_table_vec", "ing_dt_ing_calendar_month_t_desde_primera_adm_dedup_wide2_cens", "ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens", "lcs_dist_month", "lcs_dist_month_c", "lcs_dist_quarter", "lcs_dist_quarter_c", "max_clusters", "multinom_pivot_wider", "new_labels", "new_labels2", "num_cores", "om_dist_month", "om_dist_month_c", "om_dist_quarter", "om_dist_quarter_c", "oneway_anova_effect_size", "pamRange_month_lcs", "pamRange_month_lcs2", "pamRange_month_om", "pamRange_month_om2", "pamRange_quarter_lcs", "pamRange_quarter_lcs2", "pamRange_quarter_om", "pamRange_quarter_om2", "ratio_plot", "resultados_list", "seq_plot_hc_om2_m", "seq_plot_hc_om3_m", "seq_plot_pam_om2_m", "seq_plot_pam_om2_q", "seq_plot_pam_om3_q", "seq_plot_pam_om4_q", "seq_plot_pam_om6_q", "sil_hc_om_clus2_m", "sil_hc_om_clus2_q", "sil_hc_om_clus3_m", "sil_pam_om_clus2_m", "sil_pam_om_clus2_q", "sil_pam_om_clus3_q", "sil_pam_om_clus4_q", "sil_pam_om_clus6_q", "smd_bin", "States_Wide.seq_month_t_prim_adm", "States_Wide.seq_month_t_prim_adm_cens", "States_Wide.seq_month_t_prim_adm_RM_cens", "States_Wide.seq_quarter_t_prim_adm", "States_Wide.seq_quarter_t_prim_adm_cens", "States_Wide.seq_quarter_t_prim_adm_RM_cens", "seqcount_t", "seq_mean_t", "seqtrate_t")# Filtrar el entorno global y eliminar todo lo que NO esté en ese vectorrm(list =setdiff(ls(), mis_objetos)); gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 4480875 239.4 16489442 880.7 20611802 1100.8
Vcells 358216250 2733.0 2136969489 16303.8 2569386859 19602.9
Tiempo que demora esta sección: 0.1 minutos
Resultados
1. Trimestral
1.1. PAM (OM), sol 4 cluster- diagnósticos
Código
invisible("Me da buena: 0,61 en promedio. El problema está con 5710 y 6036 que son pésimos")sil_pam_om_clus4_q_nostd<-silhouette(as.integer(pamRange_quarter_om$clustering$cluster4), as.dist(dist_quarter_om))# Crear etiquetas personalizadascluster_labels <-paste0("Cluster ", seq_along(attr(summary(sil_pam_om_clus4_q_nostd)$clus.avg.widths, "dimnames")[[1]]), ":\nAWS ", sprintf("%1.2f",summary(sil_pam_om_clus4_q_nostd)$clus.avg.widths))# Graficar con etiquetas personalizadasfviz_silhouette( sil_pam_om_clus4_q_nostd, lab.clusters = cluster_labels, # Etiquetas personalizadas para los clústeresprint.summary=F) +scale_fill_grey(start =0.2, end =0.8, labels = cluster_labels) +# Escala de grisesscale_color_grey(start =0.2, end =0.8, labels = cluster_labels)+# Escala de grises para los bordeaggtitle(NULL)+labs(y="Ancho medio de la silueta", x="Conglomerados")# Elimina el título
cat("de quienes pertenecen al grupo de un semestre")
de quienes pertenecen al grupo de un semestre
Código
length(subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, clus_pam_om4=="6522, Un semestre TSM(1)")$run)
[1] 399
Código
diag_pam_om4_6522<-df_filled2 |> dplyr::filter(run %in%subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, clus_pam_om4=="6522, Un semestre TSM(1)")$run) |> dplyr::select(run, diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8, diag9, diag10, diag11, fecha_egreso_rec_fmt, estab_homo) |> dplyr::group_by(run) |> dplyr::filter(row_number() !=1) |># Elimina la primera observación de cada run dplyr::mutate(all_diags =paste(na.omit(c(diag1, diag2, diag3, diag4, diag5, diag6, diag7, diag8, diag9, diag10, diag11)), collapse =", ") ) |> dplyr::summarise(all_diags =first(all_diags),fecha_egreso_rec_fmt =first(fecha_egreso_rec_fmt),estab_homo =first(estab_homo) ) |> dplyr::ungroup() |> dplyr::pull(all_diags) |># Extraer la columna all_diags como vectorstrsplit(split =", ") |># Separar cada diagnóstico por comasunlist() # Convertir la lista en un vectorinvisible("head(arrange(data.frame(table(diag_pam_om4_6522)),-Freq),20)")invisible("Para chatgpt= estos son códigos de CIE-10, descríbeme brevemente cada uno en markdown en formato 'Cód. CIE-10 (n=Freq) - [descripción] '")cat("número de diagnósticos distintos")
## F329 (n=134) - Episodio depresivo no especificado # F322 (n=103) - Episodio depresivo grave sin síntomas psicóticos # F609 (n=96) - Otros trastornos específicos de la personalidad # F603 (n=95) - Trastorno de la personalidad emocionalmente inestable, tipo límite # F209 (n=66) - Esquizofrenia no especificada # F192 (n=53) - Trastornos mentales y del comportamiento debidos al uso de múltiples drogas y al uso de otras sustancias psicoactivas # F319 (n=49) - Trastorno depresivo recurrente, episodio actual no especificado # F200 (n=45) - Esquizofrenia paranoide # F432 (n=38) - Trastornos de adaptación # F323 (n=35) - Episodio depresivo grave con síntomas psicóticos # Z915 (n=32) - Antecedentes personales de traumatismo no clasificado en otra parte # C490 (n=29) - Neoplasia maligna del tejido conectivo y de los tejidos blandos, de localización no especificada # G909 (n=29) - Trastorno del sistema nervioso no especificado # F431 (n=28) - Reacción al estrés agudo y trastornos de adaptación # E669 (n=26) - Obesidad, no especificada# Z511 (n=24) - Atención sanitaria para radioterapia # F29X (n=21) - Psicosis no orgánica no especificada # F419 (n=20) - Trastorno de ansiedad no especificado # C901 (n=16) - Leucemia de células plasmáticas# F199 (n=16) - Trastornos mentales y de compurtamiento por el uso de múltiples drogas y uso de otras sustancias piscoactivas
Tiempo que demora esta sección: 0 minutos
Entre quienes se encuentren en un semestre en el sistema por TSM y presentan un segundo episodio, de las 20 condiciones más frecuentes, al menos el 50% se caracteriza por episodios depresivos no especificado (F329), trastornos de la personalidad tipo límite (F603), otros trastornos específicos de la personalidad (F609) y episodios depresivos graves sin síntomas psicóticos (F322) y esquizofrenia no especificada (F209).
Código
# Trastornos mentales orgánicos: F00.0-F09.9 (F000 a F099)in_organic <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^F0[0-9]", codigo)})# Trastornos por uso de sustancias: F10.0-F19.9 (F10 a F19)in_substance <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^F1[0-9]", codigo)})# Esquizofrenia: F20.0-F20.9 (F20)in_esquizofrenia <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^F20", codigo)})# Otros trastornos psicóticos no afectivos: F21.0-F29.9 (F21 a F29)in_otros_psicoticos <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^F2[1-9]", codigo)})# Trastornos bipolares: F30.0-F31.9 (F30 o F31)in_bipolares <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^(F30|F31)", codigo)})# Trastornos depresivos y otros del estado de ánimo: F32.0-F39.9 (F32 a F39)in_depresivos <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^F3[2-9]", codigo)})# Trastornos de ansiedad: F40.0-F49.9 (F40 a F49)in_ansiedad <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^F4[0-9]", codigo)})# Trastornos de la personalidad: F60.0-F69.9 (F60 a F69)in_personalidad <-Vectorize(function(codigo) {if (is.na(codigo)) return(FALSE)grepl("^F6[0-9]", codigo)})agregar_columnas_icd <-function(df) {# Especificamos las columnas de diagnóstico diag_cols <-paste0("diag", 1:11) df |>mutate(organic =rowSums(across(all_of(diag_cols), ~in_organic(.))) >0,substance =rowSums(across(all_of(diag_cols), ~in_substance(.))) >0,esquizofrenia =rowSums(across(all_of(diag_cols), ~in_esquizofrenia(.))) >0,otros_psicoticos =rowSums(across(all_of(diag_cols), ~in_otros_psicoticos(.))) >0,bipolares =rowSums(across(all_of(diag_cols), ~in_bipolares(.))) >0,depresivos =rowSums(across(all_of(diag_cols), ~in_depresivos(.))) >0,ansiedad =rowSums(across(all_of(diag_cols), ~in_ansiedad(.))) >0,personalidad =rowSums(across(all_of(diag_cols), ~in_personalidad(.))) >0 )}#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_cat("2025-03-01: general")
$chisq
Pearson's Chi-squared test
data: df
X-squared = 88.475, df = 6, p-value < 2.2e-16
$percentages
glosa_pueblo_originario 6623, Un trimestre, TSM(4) 6612, Un trimestre, TUS(3)
MAPUCHE 0.02 0.04
NINGUNO 0.98 0.96
RAPA NUI (PASCUENSE) 0.00 0.00
6522, Un semestre TSM(1) 6574, Comorbilidad un trimestre(2)
0.02 0.01
0.96 0.98
0.01 0.01
Código
invisible("No exclusivo de la clasificación por conglomerados, sino de las trayectorias en general y las matrices de distancia")#https://cran.r-project.org/web/packages/TraMineRextras/TraMineRextras.pdfseqcompare_sex_quarter_om<-seqCompare(States_Wide.seq_quarter_t_prim_adm, group=ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$glosa_sexo, method="OM", sm=costmatrix_quarter, seed=2125, s=1e1, stat="all")# LRT p-value BIC diff. Bayes Factor (Avg) Bayes Factor (From Avg BIC)# [1,] 3.723906 0.07729539 0.7281733 2.140152 1.439199seqcompare_sex_quarter_lcs<-seqCompare(States_Wide.seq_quarter_t_prim_adm, group=ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$glosa_sexo, method="LCS", sm=costmatrix_quarter, seed=2125, s=1e1, stat="all")# LRT p-value BIC diff. Bayes Factor (Avg) Bayes Factor (From Avg BIC)# [1,] 3.537305 0.08879071 0.5415731 2.120956 1.310995seqcompare_ppoo_quarter_om<-seqCompare(States_Wide.seq_quarter_t_prim_adm, group=ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$inclusivo_real_historico, method="OM", sm=costmatrix_quarter, seed=2125, s=1e1, stat="all")# Currently seqLRT supports only 2 groups!# LRT p-value BIC diff. Bayes Factor (Avg) Bayes Factor (From Avg BIC)# [1,] 3.321524 0.08714022 0.325792 1.479436 1.176914seqcompare_ppoo_quarter_lcs<-seqCompare(States_Wide.seq_quarter_t_prim_adm, group=ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$inclusivo_real_historico, method="LCS", sm=costmatrix_quarter, seed=2125, s=1e1, stat="all")# LRT p-value BIC diff. Bayes Factor (Avg) Bayes Factor (From Avg BIC)# [1,] 3.214252 0.09793892 0.21852 1.456759 1.115452# glosa_pueblo_originario_rec<- df_filled2[,c("run","glosa_pueblo_originario")] |> dplyr::left_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run", "clus_pam_om4","factor_inclusivo_real_hist_mas_autperc")], by="run", multiple="first") |> dplyr::mutate(glosa_pueblo_originario_rec= dplyr::case_when(glosa_pueblo_originario=="NINGUNO"& factor_inclusivo_real_hist_mas_autperc!="00"~"DESCONOCIDO", glosa_pueblo_originario=="NINGUNO"~"NINGUNO", grepl("RAPA|MAPU", glosa_pueblo_originario)~"RAPA NUI o MAPUCHE", T~"RESTO")) |>mutate(pueblo_bin =case_when( glosa_pueblo_originario_rec =="RAPA NUI o MAPUCHE"~"Mapuche o Rapa Nui", glosa_pueblo_originario_rec %in%c("DESCONOCIDO", "NINGUNO", "RESTO") ~"Otro o ninguno" ) ) |>pull(pueblo_bin)seqcompare_ppoo_quarter_om_ppoo_rec<-seqCompare(States_Wide.seq_quarter_t_prim_adm, group=glosa_pueblo_originario_rec, method="OM", sm=costmatrix_quarter, seed=2125, s=1e1, stat="all")
Error: ‘with.missing’ must be TRUE when ‘seqdata’ or ‘refseq’ contain missing values
Tiempo que demora esta sección: 0.1 minutos
Generamos un gráfico de PPOO por cada conglomerado.
Código
ppoo_clus_pre_pam_om4_q<- df_filled2[,c("run","glosa_pueblo_originario")] |> dplyr::left_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run", "clus_pam_om4","factor_inclusivo_real_hist_mas_autperc")], by="run", multiple="first") |> dplyr::mutate(glosa_pueblo_originario_rec= dplyr::case_when(glosa_pueblo_originario=="NINGUNO"& factor_inclusivo_real_hist_mas_autperc!="00"~"DESCONOCIDO", T~glosa_pueblo_originario)) |> janitor::tabyl(glosa_pueblo_originario_rec, clus_pam_om4) |> janitor::adorn_percentages("row")reshape2::melt(ppoo_clus_pre_pam_om4_q, id.vars ="glosa_pueblo_originario_rec") |> dplyr::mutate(glosa_pueblo_originario_rec= dplyr::recode(glosa_pueblo_originario_rec, "OTRO (ESPECIFICAR)"="OTRO(n=86)", #"OTRO(n=77)", "RAPA NUI (PASCUENSE)"="RAPA NUI(n=37)", #"RAPA NUI(n=34)", "YAGÁN (YÁMANA)"="YAGÁN(n=2)",#"YAGÁN(n=2)","AYMARA"="AYMARA(n=15)",#"AYMARA(n=13)","COLLA"="COLLA(n=6)",#"COLLA(n=6)","DIAGUITA"="DIAGUITA(n=3)",#"DIAGUITA(n=3)","KAWÉSQAR"="KAWÉSQAR(n=4)",#"KAWÉSQAR(n=4)","MAPUCHE"="MAPUCHE(n=299)",#"MAPUCHE(n=255)","DESCONOCIDO"=".DESCONOCIDO(n=2.353)",#".DESCONOCIDO(n=1.985)","NINGUNO"=".NINGUNO(n=10.425)"#".NINGUNO(n=9.156)" )) |>ggplot(aes(x = glosa_pueblo_originario_rec, y = value, fill = variable)) +geom_bar(stat ="identity", position ="fill") +scale_fill_manual(values =c("#D2B48C", "#E27A5B", "#708090", "#6B8E23")) +labs(title =NULL,x ="Grupo Étnico",y ="Proporción de reportes",fill ="Grupos") +# Cambia el título de la leyenda a "Grupos"theme_minimal() +theme(axis.text.y =element_text(size =12), # Tamaño de las etiquetas de los grupos étnicosaxis.text.x =element_text(size =12), # Tamaño de las etiquetas del eje Xaxis.title.x =element_text(size =14), # Tamaño del título del eje Xaxis.title.y =element_text(size =14), # Tamaño del título del eje Yplot.title =NULL, # Tamaño y estilo del título del gráficolegend.title =element_text(size =14, margin =margin(b =-.1)), # Tamaño del título de la leyendalegend.spacing.y =unit(1.5, "lines"),legend.box.spacing =unit(0.5, "lines"), # Controla el espacio entre la leyenda y el gráficolegend.margin =margin(5, 5, 5, 5), legend.key.height =unit(1, "cm"), legend.text =element_text(size =12) # Tamaño del texto de la leyenda ) +coord_flip() # Hacer el gráfico horizontalggsave("_figs/grafico_ancho_achatado_pam_om4_q_25.png", width =10, height =5, dpi=1000)ppoo_clus_pre_pam_om4_q_rapanui<- df_filled2[,c("run","glosa_pueblo_originario")] |> dplyr::left_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run", "clus_pam_om4","factor_inclusivo_real_hist_mas_autperc")], by="run", multiple="first") |> dplyr::mutate(glosa_pueblo_originario_rec= dplyr::case_when(glosa_pueblo_originario=="NINGUNO"& factor_inclusivo_real_hist_mas_autperc!="00"~"DESCONOCIDO", glosa_pueblo_originario=="NINGUNO"~"NINGUNO", glosa_pueblo_originario=="RAPA NUI (PASCUENSE)"~"RAPA NUI", T~"RESTO")) |> janitor::tabyl(glosa_pueblo_originario_rec, clus_pam_om4) |> janitor::adorn_percentages("row")ppoo_clus_pre_pam_om4_q_rapanuicat("origen?")df_filled2[,c("run","glosa_pueblo_originario")] |> dplyr::left_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run", "clus_pam_om4","factor_inclusivo_real_hist_mas_autperc", "codigo_region_rec_base")], by="run", multiple="first") |> dplyr::mutate(glosa_pueblo_originario_rec= dplyr::case_when(glosa_pueblo_originario=="NINGUNO"& factor_inclusivo_real_hist_mas_autperc!="00"~"DESCONOCIDO", glosa_pueblo_originario=="NINGUNO"~"NINGUNO", glosa_pueblo_originario=="RAPA NUI (PASCUENSE)"~"RN", T~"RESTO")) |> dplyr::filter(glosa_pueblo_originario_rec=="RN") |> janitor::tabyl(codigo_region_rec_base)
glosa_pueblo_originario_rec 6623, Un trimestre, TSM(4)
DESCONOCIDO 0.7403315
NINGUNO 0.7310312
RAPA NUI 0.1891892
RESTO 0.7301205
6612, Un trimestre, TUS(3) 6522, Un semestre TSM(1)
0.13047174 0.09562261
0.09870504 0.12853717
0.13513514 0.54054054
0.12289157 0.13012048
6574, Comorbilidad un trimestre(2)
0.03357416
0.04172662
0.13513514
0.01686747
origen? codigo_region_rec_base n percent
RM 22 0.5945946
noRM 15 0.4054054
Eventos hospitalarios y reportes de pertenencia a PPOO MINSAL, por conglomerados
Tiempo que demora esta sección: 0 minutos
1.1.1. Trayectorias
Vemos los gráficos de las trayectorias
Código
categories_pam_om4_q<-attr(States_Wide.seq_quarter_t_prim_adm_cens, "labels")new_labels <- categories_pam_om4_qnew_labels[which(categories_pam_om4_q =="Otras causas")] <-"Otras\ncausas"#new_labels[which(categories == "Consumo\nde sustancias")] <- "Consumo de\nsustancias"# Creamos un vector con las columnas llenando con NA si faltan valoressil_pam_om_clus4_q <-wcSilhouetteObs(as.dist(dist_quarter_om), pamRange_quarter_om$clustering$cluster4, measure="ASW")seq_plot_pam_om4_q <-ggseqiplot(States_Wide.seq_quarter_t_prim_adm_cens, group= ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4,facet_ncol=2, facet_nrow=2, sortv=sil_pam_om_clus4_q) +theme(legend.position ="none")+labs(x="Trimestres", y="# IDs de usuarios")+#guides(fill = guide_legend(nrow = 1))+theme(panel.spacing =unit(0.1, "lines"), # Reduce el espaciado entre los panelesaxis.text.y =element_text(size =15), # Tamaño de las etiquetas de los grupos étnicosaxis.text.x =element_text(size =15), # Tamaño de las etiquetas del eje Xaxis.title.x =element_text(size =15), # Tamaño del título del eje Xaxis.title.y =element_text(size =15, margin =margin(r =-10)),#,margin = margin(l = -10)),strip.text =element_text(size =11, margin =margin(b =-15)),legend.text =element_text(size =15),legend.spacing.x =unit(0.1, 'cm'), # Alinea el título de la leyenda hacia la izquierdalegend.box.margin =margin(t =0, r =0, b =0, l =-50),legend.position ="bottom", legend.justification ="left",panel.spacing.y =unit(0.5, "lines"),strip.placement ="outside", # Para colocar las tiras fuera de los ejesstrip.background =element_blank() # Elimina el fondo para que parezca más espacioso#legend.key.size = unit(1.5, "lines"), # Aumenta el tamaño de los símbolos en la leyenda )+guides(fill =guide_legend(nrow =1)) +scale_fill_manual(labels = new_labels, values=c("#E2725B", "#556B2F", "#D2B48C",#"#8B4513","#FFFFFF","#808080","#000000"))+scale_color_manual(labels = new_labels, values=c("#E2725B", "#556B2F", "#D2B48C",#"#8B4513","#FFFFFF","#808080","#000000"))seq_plot_pam_om4_q ggsave(filename="_figs/clusters_pam_om4_q_mod_25.png", seq_plot_pam_om4_q, width =9.5, height =5.5, dpi=1000)
Trayectorias de hospitalización, orden de sujetos según el primer estado observado y su duración, representando a cada individuo como una línea en el gráfico (observaciones ordenadas de acuerdo a ASW)
Tiempo que demora esta sección: 0.3 minutos
Código
seq_plot2_pam_om4_q <-ggseqdplot(States_Wide.seq_quarter_t_prim_adm_cens, group= ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4,facet_ncol=2, facet_nrow=2) +theme(legend.position ="none")+# Colocar la leyenda abajolabs(x="Trimestres", y="Frecuencia relativa de estados")+theme(panel.spacing =unit(0.1, "lines"),axis.text.y =element_text(size =15), # Tamaño de las etiquetas de los grupos étnicosaxis.text.x =element_text(size =15), # Tamaño de las etiquetas del eje Xaxis.title.x =element_text(size =15), # Tamaño del título del eje Xaxis.title.y =element_text(size =15, margin =margin(r =-5)),strip.text =element_text(size =11),panel.spacing.y =unit(0.5, "lines"),strip.placement ="outside", # Para colocar las tiras fuera de los ejesstrip.background =element_blank() # Elimina el fondo para que parezca más espacioso#legend.key.size = unit(1.5, "lines"), # Aumenta el tamaño de los símbolos en la leyenda ) # Colocar la leyenda abajoseq_plot2_pam_om4_qggsave("_figs/clusterspam_om42_q_mod_25.png",seq_plot2_pam_om4_q, width =8.5, height =5.5, dpi=1000)table_data_pam_om4_q <-sprintf("%1.2f",pamRange_quarter_om$stats[3,])table_data_pam_om4_q <-as.data.frame(t(table_data_pam_om4_q))colnames(table_data_pam_om4_q)<-attr(pamRange_quarter_om$stats, "name")table_data_pam_om4_q |> knitr::kable()
PBC
HG
HGSD
ASW
ASWw
CH
R2
CHsq
R2sq
HC
0.49
0.66
0.65
0.58
0.59
948.53
0.30
1312.64
0.37
0.20
Trayectorias de hospitalización, frecuencia relativa de estados en un gráfico de barras apiladas por trimestre.
Trayectorias de hospitalización, frecuencia relativa de estados en un gráfico de barras apiladas por trimestre.
Tiempo que demora esta sección: 0.1 minutos
De este modo, presenta el cambio agregado en la distribución de estados a lo largo del tiempo, sin considerar las secuencias individuales.
Código
cat("Definimos las observaciones que tienen siluetas negativas")sil_neg_pam_om_clus4_q <-which(sil_pam_om_clus4_q<0)cat("A qué conglomerados pertenecen?")table(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[sil_neg_pam_om_clus4_q, "clus_pam_om4"])ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$rn<-1:nrow(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens)
Definimos las observaciones que tienen siluetas negativasA qué conglomerados pertenecen?clus_pam_om4
6623, Un trimestre, TSM(4) 6612, Un trimestre, TUS(3)
2 0
6522, Un semestre TSM(1) 6574, Comorbilidad un trimestre(2)
151 0
Tiempo que demora esta sección: 0 minutos
1.1.2.Exploración transiciones
1.1.2.a Transiciones- RM y no RM
Tasas de transición no RM a RM y viceversa
Código
invisible("Tasas de transición no RM a RM y viceversa")trim_tasa_pam_om4_q_cens_cnt<-seqcount_t(States_Wide.seq_quarter_t_prim_adm_RM_cens, group=ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4) |> dplyr::filter(count>0) |> dplyr::mutate(trans =paste0(from,"_", to)) |> dplyr::mutate(across(c("from","to"),~gsub("\\[->\\s*|\\s*->\\s*\\]|\\[|\\]", "", .))) trim_tasa_pam_om4_q_cens_rate<-seqtrate_t(States_Wide.seq_quarter_t_prim_adm_RM_cens, group=ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4) |> dplyr::filter(rate>0) |> dplyr::mutate(trans =paste0(from,"_", to)) |> dplyr::mutate(across(c("from","to"),~gsub("\\[->\\s*|\\s*->\\s*\\]|\\[|\\]", "", .)))
Tiempo que demora esta sección: 0 minutos
Código
trim_tasa_pam_om4_q_cens_rate |> dplyr::left_join(trim_tasa_pam_om4_q_cens_cnt, by=c("from"="from", "glosa_sexo"="glosa_sexo","to"="to")) |> dplyr::rename("recuento"="count") |> dplyr::filter(from %in%c("RM", "noRM")) |>ggplot(aes(x = from, y = to, fill = rate, size=log(recuento+1))) +geom_tile() +coord_flip()+scale_fill_gradient(low ="white", high ="blue") +# Ajusta la escala de colores según tus preferenciaslabs(title ="Tasas de transición, Trimestre (s/censura)",x ="Desde",y ="Hacia",fill ="Rate") +theme_minimal() +facet_wrap(~glosa_sexo)+theme(axis.text.x =element_text(angle =45, hjust =1))+geom_text(aes(label =sprintf("%1.2f", rate), size =log(recuento+1)*.5), color ="black")invisible("Hay muy pocos casos que se entrecruzan entre noRM y RM (fuera de la diagnonal)")
Porcentajes de transición no-RM y RM por cada cluster
Tiempo que demora esta sección: 0.1 minutos
Hay muy pocos casos que se entrecruzan entre noRM y RM (fuera de la diagnonal)
trim_tasa2_pam_om4_q_cens_rate |> dplyr::left_join(trim_tasa2_pam_om4_q_cens_cnt, by=c("from"="from", "glosa_sexo"="glosa_sexo","to"="to")) |> dplyr::rename("recuento"="count") |>#dplyr::filter(from %in% c("RM", "noRM")) |> ggplot(aes(x = from, y = to, fill = rate, size=log(recuento+1))) +geom_tile() +coord_flip()+scale_fill_gradient(low ="white", high ="blue") +# Ajusta la escala de colores según tus preferenciaslabs(title ="Tasas de transición, Trimestre (s/censura)",x ="Desde",y ="Hacia",fill ="Rate") +theme_minimal() +facet_wrap(~glosa_sexo)+theme(axis.text.x =element_text(angle =45, hjust =1))+geom_text(aes(label =sprintf("%1.2f", rate), size =log(recuento+1)*.5), color ="black")
Porcentajes de transición, transiciones posteriores, por cada cluster
Tiempo que demora esta sección: 0.1 minutos
1.1.2.c Tiempo promedio por cluster
Código
seq_mean_t(States_Wide.seq_quarter_t_prim_adm_cens, group=ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4) |> data.table::as.data.table(keep.rowname=T) |> dplyr::mutate(rn=gsub("\\d", "", rn)) |>ggplot(aes(x=rn, fill= factor_inclusivo, y=Mean))+geom_bar(width =1, stat ="identity") +theme_minimal() +facet_wrap(~factor_inclusivo)+labs(title =NULL,x =NULL,y =NULL) +scale_fill_manual(values =c("#70809090", "#6B8E2380", "#E27A5B","#D2B48C")) +coord_flip()+theme(#axis.text.x = element_blank(),#axis.text.y = element_blank(),panel.grid =element_blank()) +# scale_fill_brewer(palette = "Pastel1", labels=c("Sin\nautoidentificación\nni reconocimiento", "Autoidentificación\nsin reconocimiento", "Ambas")) +geom_text(aes(label =round(Mean,1)), position =position_stack(vjust =0.5), size =3.5, # Ajusta el tamaño de la fuente aquícolor ="black", # Color del textofamily ="sans", # Puedes cambiar la fuente si lo deseasbackground =element_rect(fill ="white", color =NA)) +# Fondo blancotheme(legend.title =element_blank())invisible("No me aporta mucho")
Tiempo promedio en cada estado por estatus PPOO (Trimestral c/censura)
Tiempo que demora esta sección: 0.1 minutos
Observamos que aquellos en el conglomerado que se encuentra un trimestre, la duración de los ingresos relacionados con trastornos de salud mental son en promedio más largos.
1.1.3. Comparación variables
1.1.3.a. Comparación covariables- PPOO
Código
ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens|> dplyr::count(clus_pam_om4, factor_inclusivo_real_hist_mas_autperc)|> dplyr::group_by(clus_pam_om4)|> dplyr::mutate(n_prop =paste0(n, " (",scales::percent(n /sum(n), accuracy=.1),")"))|> dplyr::select(-n)|> tidyr::pivot_wider(names_from = clus_pam_om4, values_from = n_prop, values_fill ="0")|> dplyr::mutate(factor_inclusivo_real_hist_mas_autperc =factor(factor_inclusivo_real_hist_mas_autperc, levels =c("11", "10", "00"), labels=c("Se identifica/hay reconocimeinto", "No se identifica/hay reconocimiento", "No se identifica/no pertenece")))|> dplyr::arrange(factor_inclusivo_real_hist_mas_autperc)|> (\(df) {if (interactive()) {mutate(df, across(-factor_inclusivo_real_hist_mas_autperc, ~gsub("%", "", gsub("\\.", ",", .))))|> rio::export("clipboard")} knitr::kable(df, caption="Porcentajes por fila, conglomerado vs. Pertenencia/identificación + Reconocimento CONADI PPOO") })()
Porcentajes por fila, conglomerado vs. Pertenencia/identificación + Reconocimento CONADI PPOO
factor_inclusivo_real_hist_mas_autperc
6623, Un trimestre, TSM(4)
6612, Un trimestre, TUS(3)
6522, Un semestre TSM(1)
6574, Comorbilidad un trimestre(2)
Se identifica/hay reconocimeinto
443 (8.4%)
75 (10.0%)
27 (6.8%)
17 (7.4%)
No se identifica/hay reconocimiento
596 (11.4%)
92 (12.3%)
45 (11.3%)
21 (9.2%)
No se identifica/no pertenece
4209 (80.2%)
583 (77.7%)
327 (82.0%)
191 (83.4%)
Tiempo que demora esta sección: 0 minutos
Vemos las categorías de clasificación de PPOO según autopercepción (en MINSAL y en RSH) y reconocimiento CONADI.
Porcentajes por fila, conglomerado vs. Pertenencia/identificación + Reconocimento CONADI PPOO
Conglomerados
No se identifica/no pertenece
Hay reconocimeinto
6623, Un trimestre, TSM(4)
4209 (80.2%)
1039 (19.8%)
6612, Un trimestre, TUS(3)
583 (77.7%)
167 (22.3%)
6522, Un semestre TSM(1)
327 (82.0%)
72 (18.0%)
6574, Comorbilidad un trimestre(2)
191 (83.4%)
38 (16.6%)
Tiempo que demora esta sección: 0 minutos
1.1.3.b. Comparación covariables- Mortalidad
Código
# invisible("No hay nada, el tiempo promedio de censura es similar")ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens |> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)) |> janitor::tabyl(clus_pam_om6,death_time_rec) |> dplyr::mutate(`1`=paste0(`1`," (", scales::percent(`1`/(`0`+`1`), accuracy=.1),")")) |> dplyr::left_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens |> dplyr::group_by(clus_pam_om6) |> dplyr::summarise(mean=sprintf("%1.1f",mean(cens_time))), by="clus_pam_om6") |> dplyr::select(-`0`)|> (\(df) {if (interactive()) {mutate(df, across(-1, ~gsub("%", "", gsub("\\.", ",", .))))|>select(2)|>t()|> rio::export("clipboard")} knitr::kable(df, "markdown", col.names=c("Conglomerado","Mortalidad observada", "Promedio"), caption="Post-hoc, conglomerado vs. Mortalidad y tiempo a censura") })()
Post-hoc, conglomerado vs. Mortalidad y tiempo a censura
Conglomerado
Mortalidad observada
Promedio
6623, Un trimestre, TSM(5)
58 (1.2%)
17.9
6612, Un trimestre, TUS(4)
19 (2.5%)
18.2
6522, Un semestre TSM(2)
8 (2.2%)
17.9
6624, TSM, 1 año después, otras causas(6)
5 (2.1%)
17.8
6574, Comorbilidad un trimestre(3)
6 (2.7%)
18.1
6268, TSM, 1 año después, TSM(1)
3 (1.7%)
18.1
Tiempo que demora esta sección: 0 minutos
Código
ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens |> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)) |> janitor::tabyl(death_time_rec,clus_pam_om4) |> janitor::chisq.test(correct=T)#X-squared = 11.377, df = 3, p-value = 0.009854chisq_cramerv(with(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens|> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)), table(death_time_rec , clus_pam_om4)))fisher.test(with(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens|> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)), table(death_time_rec , clus_pam_om4), simulate.p.value=T, B=1e5))#p-value = 0.008208message("Descartando valores negativos en sil width")chisq_cramerv(with(subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, !rn %in% sil_neg_pam_om_clus4_q)|> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)), table(death_time_rec , clus_pam_om4)))# $chisq_statistic# [1] "11.31"# # $chisq_df# df # 3 # # $chisq_p_value# [1] "0.0101"# # $cramers_v# [1] "0.04"fisher.test(with(subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, !rn %in% sil_neg_pam_om_clus4_q)|> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)), table(death_time_rec , clus_pam_om4), simulate.p.value=T, B=1e5))#p-value = 0.007751invisible("no se basa en la distribución chi-cuadrado. Fisher se basa en permutaciones exactas, por lo que no se calculan df.")#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_##_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_tab_cl_mortalidad_pam_om4_q<- ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens |> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)) |> janitor::tabyl(death_time_rec,clus_pam_om4) |>as.matrix(ncol=2)labels_pam_om4_q <-c("6623, Un trimestre, TSM(4)","6612, Un trimestre, TUS(3)","6522, Un semestre TSM(1)","6574, Comorbilidad un trimestre(2)")chisq.posthoc.test(with(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens|> dplyr::mutate(death_time_rec=ifelse(death_time==20,0,1)), table(death_time_rec , clus_pam_om4), simulate.p.value=T, B=1e5))# Realizar el análisis y crear la tabla directamentepairwise.prop.test(t(tab_cl_mortalidad_pam_om4_q[,2:5]), p.adjust.method ="holm")$p.value|>as.table() |>as.data.frame() |>rename(Grupo_1 = Var1, Grupo_2 = Var2, p_value = Freq) |>filter(!is.na(p_value)) |>mutate(Grupo_1 = labels_pam_om4_q[as.numeric(Grupo_1)],Grupo_2 = labels_pam_om4_q[as.numeric(Grupo_2)],p_value =ifelse(p_value <.001, "<.001", sprintf("%1.3f",p_value)) ) |>kable(col.names =c("Grupo 1", "Grupo 2", "Valor p ajustado"),align ="l",caption="Corrección parcial por comparaciones múltiples (Holm–Bonferroni)" )
Pearson's Chi-squared test
data: janitor::tabyl(dplyr::mutate(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, death_time_rec = ifelse(death_time == 20, 0, 1)), death_time_rec, clus_pam_om4)
X-squared = 11.377, df = 3, p-value = 0.009854
$chisq_statistic
[1] "11.38"
$chisq_df
df
3
$chisq_p_value
[1] "0.0099"
$cramers_v
[1] "0.04"
Fisher's Exact Test for Count Data
data: with(dplyr::mutate(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, death_time_rec = ifelse(death_time == 20, 0, 1)), table(death_time_rec, clus_pam_om4), simulate.p.value = T, B = 1e+05)
p-value = 0.008208
alternative hypothesis: two.sided
$chisq_statistic
[1] "11.31"
$chisq_df
df
3
$chisq_p_value
[1] "0.0101"
$cramers_v
[1] "0.04"
Fisher's Exact Test for Count Data
data: with(dplyr::mutate(subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, !rn %in% sil_neg_pam_om_clus4_q), death_time_rec = ifelse(death_time == 20, 0, 1)), table(death_time_rec, clus_pam_om4), simulate.p.value = T, B = 1e+05)
p-value = 0.007751
alternative hypothesis: two.sided
Dimension Value 6623, Un trimestre, TSM(4) 6612, Un trimestre, TUS(3)
1 0 Residuals 3.346154 -2.491149
2 0 p values 0.006555 0.101865
3 1 Residuals -3.346154 2.491149
4 1 p values 0.006555 0.101865
6522, Un semestre TSM(1) 6574, Comorbilidad un trimestre(2)
1 -1.293403 -1.429422
2 1.000000 1.000000
3 1.293403 1.429422
4 1.000000 1.000000
Corrección parcial por comparaciones múltiples (Holm–Bonferroni)
Grupo 1
Grupo 2
Valor p ajustado
6623, Un trimestre, TSM(4)
6623, Un trimestre, TSM(4)
0.047
6612, Un trimestre, TUS(3)
6623, Un trimestre, TSM(4)
0.654
6522, Un semestre TSM(1)
6623, Un trimestre, TSM(4)
0.654
6612, Un trimestre, TUS(3)
6612, Un trimestre, TUS(3)
1.000
6522, Un semestre TSM(1)
6612, Un trimestre, TUS(3)
1.000
6522, Un semestre TSM(1)
6522, Un semestre TSM(1)
1.000
Tiempo que demora esta sección: 0 minutos
Código
# Cargar las librerías necesariaslibrary(survival)library(ggplot2)# Crear la variable de supervivenciasurv_obj_4c <-Surv(time= ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$death_time,event =ifelse(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$death_time==20,0,1))# Realizar el análisis de Log-Rank (survdiff)surv_diff_4c <-survdiff(surv_obj_4c ~ clus_pam_om4,data = ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens)cat("Diferencia multiplicativa entre cluster TSM 1 trimestre y el resto")
Diferencia multiplicativa entre cluster TSM 1 trimestre y el resto
Pairwise comparisons using Log-Rank test
data: subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, and clus_pam_om4 !rn %in% sil_neg_pam_om_clus4_q) and clus_pam_om4
6623, Un trimestre, TSM(4)
6612, Un trimestre, TUS(3) 0.028
6522, Un semestre TSM(1) 0.438
6574, Comorbilidad un trimestre(2) 0.348
6612, Un trimestre, TUS(3)
6612, Un trimestre, TUS(3) -
6522, Un semestre TSM(1) 1.000
6574, Comorbilidad un trimestre(2) 1.000
6522, Un semestre TSM(1)
6612, Un trimestre, TUS(3) -
6522, Un semestre TSM(1) -
6574, Comorbilidad un trimestre(2) 1.000
P value adjustment method: holm
# Crear el gráfico de Kaplan-Meier con ggplot2ggplot(km_data_4c, aes(x = time, y = surv, color = strata)) +geom_step(size =1.2) +# Curvas de supervivencia#geom_ribbon(aes(ymin = lower, ymax = upper, fill = strata), alpha = 0.2, color = NA) + # Intervalos de confianzalabs(x ="Tiempo (meses)",y ="Probabilidad de Supervivencia",color ="Grupo",fill ="Grupo" ) +theme_minimal() +ylim(c(0.9,1))+theme(legend.position ="bottom")
pairwise_chisq_gof_test(table(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$clus_pam_om4,ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens$codigo_region_rec_base), p.adjust.method="holm")|> knitr::kable("markdown", caption="Dependencia categórica sol. 4 conglomerados, por pares de categorías en RM")
Dependencia categórica sol. 4 conglomerados, por pares de categorías en RM
tab_cluster_region_pam_om4_q<-ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens |> dplyr::inner_join(data_long_establecimiento_2024_std[,c("ESTAB_HOMO", "codigo_region", "nivel_de_atencion", "nivel_de_complejidad")], by =c("estab_homo_base"="ESTAB_HOMO"), multiple ="first") |> janitor::tabyl(codigo_region, clus_pam_om4) |> janitor::adorn_percentages("col") |> janitor::adorn_rounding(digits =2)#colnames(tab_cluster_region_pam_om4_q)<- c("reg", "c1", "c4", "c3", "c5", "c6", "c7", "c8", "c9", "c2")cod_reg_homo_pam_om4_q<-data.frame(codigo_region =1:16,nombre_region =c("Región de Tarapacá","Región de Antofagasta","Región de Atacama","Región de Coquimbo","Región de Valparaíso","Región del Libertador General Bernardo O'Higgins","Región del Maule","Región del Biobío","Región de La Araucanía","Región de Los Lagos","Región de Aysén del General Carlos Ibáñez del Campo","Región de Magallanes y de la Antártica Chilena","Región Metropolitana de Santiago","Región de Los Ríos","Región de Arica y Parinacota","Región de Ñuble" ),stringsAsFactors =FALSE)dplyr::mutate(tab_cluster_region_pam_om4_q, promedio_fila =rowMeans(across(2:length(colnames(tab_cluster_region_pam_om4_q))))) |> dplyr::arrange(desc(promedio_fila)) |> dplyr::left_join(cod_reg_homo_pam_om4_q, by="codigo_region") |> dplyr::select(codigo_region, nombre_region, everything()) |> dplyr::select(-promedio_fila) |> dplyr::mutate_at(3:(length(colnames(tab_cluster_region_pam_om4_q))+1),~scales::percent(.))|> knitr::kable(caption="Porcentaje por región")
Porcentaje por región
codigo_region
nombre_region
6623, Un trimestre, TSM(4)
6612, Un trimestre, TUS(3)
6522, Un semestre TSM(1)
6574, Comorbilidad un trimestre(2)
13
Región Metropolitana de Santiago
45%
35%
46%
56%
5
Región de Valparaíso
8%
13%
11%
7%
10
Región de Los Lagos
6%
19%
4%
10%
8
Región del Biobío
10%
8%
9%
9%
9
Región de La Araucanía
5%
4%
5%
5%
6
Región del Libertador General Bernardo O’Higgins
4%
3%
5%
2%
7
Región del Maule
4%
3%
3%
3%
1
Región de Tarapacá
2%
2%
2%
2%
2
Región de Antofagasta
2%
1%
4%
1%
11
Región de Aysén del General Carlos Ibáñez del Campo
pairwise_chisq_gof_test(dplyr::filter(tab_clus_macrozona_pam_om4_q,macrozona!="RM")[-1], p.adjust.method="holm")|> knitr::kable("markdown", caption="Dependencia categórica sol. 4 conglomerados, por pares de categorías en Macrozona (corrección Holm-Bonferroni)")#Groups sharing a letter are not significantlt different (alpha = 0.05)
Dependencia categórica sol. 4 conglomerados, por pares de categorías en Macrozona (corrección Holm-Bonferroni)
pairwise_chisq_gof_test(tab_clus_sexo_pam_om4_q[-1], p.adjust.method="holm")|> dplyr::mutate(p.adj= dplyr::case_when(p.adj<0.001~"<0.001",T~sprintf("%1.3f",p.adj)))|> knitr::kable("markdown", caption="Dependencia categórica sol. 4 conglomerados, por pares de categorías en Sexo (corrección Holm-Bonferroni)")
Dependencia categórica sol. 4 conglomerados, por pares de categorías en Sexo (corrección Holm-Bonferroni)
Descriptivos, edad minima de ingreso por conglomerado
clus_pam_om4
sum
6623, Un trimestre, TSM(4)
20.58 ± 4.32 ; 20 [17, 24]
6612, Un trimestre, TUS(3)
22.41 ± 4.30 ; 23 [19, 26]
6522, Un semestre TSM(1)
19.88 ± 4.28 ; 19 [16, 23]
6574, Comorbilidad un trimestre(2)
21.78 ± 4.21 ; 21 [18, 25]
Tiempo que demora esta sección: 0 minutos
Código
dt_ing_calendar_quarter_t_desde_primera_adm_dedup|>filter(quarter ==0)|>inner_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[, c("run", "clus_pam_om4")],by ="run")|> dplyr::mutate(clus_pam_om4=str_wrap(clus_pam_om4, width =20))|>group_by(clus_pam_om4)|> dplyr::summarise(mean_edad =mean(min_edad_anos),sd_edad =sd(min_edad_anos),n =n())|> dplyr::mutate(sem = sd_edad /sqrt(n),error_lower = mean_edad - sem,error_upper = mean_edad + sem)|>ggplot(aes(x = clus_pam_om4, y = mean_edad)) +geom_point(color ="black", size =3) +geom_errorbar(aes(ymin = error_lower, ymax = error_upper), width =0.2, color ="black") +labs(x ="Conglomerado (k = 4)", y ="Edad al primer evento 2018 promedio") +theme_minimal() +coord_flip() +theme(axis.text.y =element_text(size =17* .7, face ="bold"),axis.text.x =element_text(size =17* .7, face ="bold", lineheight =0.9),axis.title.x =element_text(size =16* .7, face ="bold"),axis.title.y =element_text(size =16* .7, face ="bold"),plot.margin =unit(c(0.5, 0.5, 0.5, 0.5), "cm"))ggsave("_figs/edad_minima_por_cluster_pam_om4_q_25.png", width=8, height=5, dpi=600)
Edad promedio primer ingreso con errores estándar por conglomerado
Tiempo que demora esta sección: 0.1 minutos
Código
invisible("Prueba de Levene par igualdad de varianzas")with(dt_ing_calendar_quarter_t_desde_primera_adm_dedup |> dplyr::filter(quarter ==0) |> dplyr::inner_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run","clus_pam_om4")], by="run"), car::leveneTest(min_edad_anos, clus_pam_om4))#NSanova_clus_pam_om4_q <-oneway.test(min_edad_anos ~ clus_pam_om4, data = dt_ing_calendar_quarter_t_desde_primera_adm_dedup |> dplyr::filter(quarter ==0) |> dplyr::inner_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run","clus_pam_om4")], by="run"),var.equal = T)with(dt_ing_calendar_quarter_t_desde_primera_adm_dedup |> dplyr::filter(quarter ==0) |> dplyr::inner_join(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens[,c("run","clus_pam_om4")], by="run"), oneway_anova_effect_size (min_edad_anos, clus_pam_om4))# Ver los resultados del ANOVAprint(anova_clus_pam_om4_q)#F = 49.148, num df = 3, denom df = 6622, p-value < 2.2e-16#$eta_squared#[1] 0.0217808message("Descartando valores negativos en sil width")with(dt_ing_calendar_quarter_t_desde_primera_adm_dedup |> dplyr::filter(quarter ==0) |> dplyr::inner_join(subset(ing_dt_ing_calendar_quarter_t_desde_primera_adm_dedup_wide2_cens, !rn %in% sil_neg_pam_om_clus4_q)[,c("run","clus_pam_om4")], by="run"), oneway.test(min_edad_anos ~ clus_pam_om4,var.equal = T) )#F = 48.252, num df = 3, denom df = 6469, p-value < 2.2e-16
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 3 0.4775 0.698
6622
$anova_summary
Df Sum Sq Mean Sq F value Pr(>F)
group 3 2746 915.4 49.15 <2e-16 ***
Residuals 6622 123343 18.6
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
$eta_squared
[1] 0.0217808
One-way analysis of means
data: min_edad_anos and clus_pam_om4
F = 49.148, num df = 3, denom df = 6622, p-value < 2.2e-16
One-way analysis of means
data: min_edad_anos and clus_pam_om4
F = 48.252, num df = 3, denom df = 6469, p-value < 2.2e-16
Pearson's Chi-squared test
data: df
X-squared = 30.37, df = 12, p-value = 0.002456
Fisher's Exact Test for Count Data with simulated p-value (based on
1e+05 replicates)
data: df
p-value = 0.00131
alternative hypothesis: two.sided
Porcentajes por columna, conglomerado vs. Previsión
Dependencia categórica sol. 4 conglomerados, por pares de categorías en Niv. Complejidad (corrección Holm-Bonferroni)
n
group1
group2
statistic
p
df
p.adj
p.adj.signif
5998
6623, Un trimestre, TSM(4)
6612, Un trimestre, TUS(3)
41.99007
0.00e+00
4
1.00e-07
****
5647
6623, Un trimestre, TSM(4)
6522, Un semestre TSM(1)
22.93875
1.30e-04
4
3.90e-04
***
5477
6623, Un trimestre, TSM(4)
6574, Comorbilidad un trimestre(2)
17.83641
1.33e-03
4
1.96e-03
**
1149
6612, Un trimestre, TUS(3)
6522, Un semestre TSM(1)
66.00124
0.00e+00
4
0.00e+00
****
979
6612, Un trimestre, TUS(3)
6574, Comorbilidad un trimestre(2)
18.51223
9.80e-04
4
1.96e-03
**
628
6522, Un semestre TSM(1)
6574, Comorbilidad un trimestre(2)
34.79469
5.00e-07
4
2.00e-06
****
Tiempo que demora esta sección: 0 minutos
1.1.4. Compilación comparación covariables
Código
# Definir los datos correctamentedata_pam_om4_q <-cbind.data.frame(Grupo=c("6035,\nUn trimestre,\nTSM(4)", "6025,\nUn trimestre,\nTUS(3)", "5939,\nUn semestre\nTSM(1)", "5989,\nComorbilidad un\ntrimestre(2)"), PPOO_bin =c(NA, NA, NA, NA), PPOO_sinautoid =c(NA, NA, NA, NA), PPOO_conautoid =c(NA, NA, NA, NA), Mortalidad =c(NA, NA, NA, NA), RM =c(NA, "-", NA, "+"), `Macrozona-Austral`=c(NA, NA, NA, NA), `Macrozona-Centro`=c(NA, NA, NA, NA), `Macrozona-Centro Sur`=c(NA, NA, NA, NA), `Macrozona-Norte`=c("+","-", NA, NA), `Macrozona-Sur`=c(NA, "+", NA, NA), Sexo_mujeres =c(NA, "-", NA, "-"), `Edad ingreso`=c("-", "+", "-", "+"), `Previsión-FFAA`=c("+", "-", NA, NA), `Previsión-FONASA A`=c(NA, NA, NA, NA), `Previsión-FONASA BC`=c(NA, NA, NA, NA), `Previsión-FONASA D`=c(NA, NA, NA, NA), `Previsión-ISAPRE`=c(NA, NA, NA, NA), `NivComp-Baja`=c(NA, NA, NA, NA), `NivComp-Media`=c(NA, "-", NA, NA), `NivComp-Alta`=c(NA, "+", "-", "+"))## Asegurar que los nombres de las columnas sean válidos y no haya espacios en blanco# Derretir el dataframe para que sea adecuado para ggplot2data_melt_pam_om4_q <- reshape2::melt(data_pam_om4_q, id.vars ='Grupo', variable.name ='Variable', value.name ='Asociación')# Reemplazar los NA por un valor vacíodata_melt_pam_om4_q$Asociación[is.na(data_melt_pam_om4_q$Asociación)] <-"\n"# Crear el gráfico con ggplotplot_data_melt_pam_om4_q<-data_melt_pam_om4_q |> dplyr::mutate(Variable =gsub("_", " ", Variable)) |>ggplot(aes(x = Variable, y = Grupo, fill = Asociación)) +geom_tile(color ="white", size =0.8) +scale_fill_manual(values =c("+"="#556B2F", "-"="#E2725B", "\n"="white")) +labs(title =NULL, x ="Variables", y ="Conglomerado") +theme_minimal() +theme(#axis.text.x = element_text(angle = 45, hjust = 1),panel.grid =element_blank())+theme(axis.text.y =element_text(size =17*.7, face ="bold"),#,margin = margin(l = 7)), # Tamaño de las etiquetas de los grupos étnicosaxis.text.x =element_text(size =17*.7, face ="bold"), # Tamaño de las etiquetas del eje Xaxis.title.x =element_text(size =16*.7, face ="bold"),#,margin = margin(t = -15)), # Tamaño del título del eje Xaxis.title.y =element_text(size =16*.7, face ="bold"), # Tamaño del título del eje Yplot.title =NULL, # Tamaño y estilo del título del gráficolegend.title =element_text(size =17*.7, face ="bold"), # Tamaño del título de la leyendalegend.spacing.y =unit(1.2, "lines"),legend.box.spacing =unit(0.5, "lines"), # Controla el espacio entre la leyenda y el gráficolegend.margin =margin(5, 5, 5, 5), legend.key.height =unit(1, "cm"), legend.text =element_text(size =15*.7, face ="bold") # Tamaño del texto de la leyenda ) +coord_flip()plot_data_melt_pam_om4_qggsave(filename="_figs/asociaciones_pam_om4_q_25.png", plot_data_melt_pam_om4_q, width=11*.8, height=6*.8, dpi=600)
